//============================================================================
// Name        : trjSelf.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <iterator>
#include <iomanip>
#include "AtomsXYZ.h"
#include "FstreamC.h"
#include "fxtc.h"

using namespace std;

// profile variables
extern"C" {
void openFiles(char * C_MyPDB, char * C_MyXTC, int & nato)
{
	string MyPDB(C_MyPDB);
	string MyXTC=C_MyXTC;
	ifstream ftest;
	MyPDB.erase(remove_if(MyPDB.begin(),MyPDB.end(),::isspace),MyPDB.end());
	MyXTC.erase(remove_if(MyXTC.begin(),MyXTC.end(),::isspace),MyXTC.end());
	try{
		filepdb=MyPDB;
		fpdb.open(filepdb.c_str(),ios::in);
		if(!fpdb) throw string("\n Cannot open " + filepdb + "!!\n");
		filein=MyXTC;
		ftest.open(filein.c_str(),ios::in);
		if(!ftest) throw string("\n Cannot open " + filein + "!!\n");
		ftest.close();
		finx=new FstreamC(filein,"rb");


	}catch(const string & s){
		cout << s << endl;
		exit(1);
	}

	atm=new AtomsXYZ;



	Fstream & fin=*finx;


	// read pdb file to construct topology
	vector<string> data;
	for(string str;getline(fpdb,str);){
		data.push_back(str);
	}
	// Define radii to use


	topPDB(data);
	MyTop(topPDB,false);
//	atm->setrd(MyTop);

	typedef map<string,vector<vector<int> > > mappa;
	mappa & Def=MyTop.getDef();
	vector<string> CSelect;

	int natoms=MyTop.Size();
	int nresid=MyTop.ResSize();


	atm->setDim(natoms);
	HeaderTrj header;
	nato=natoms;
// Read header of dcd file
	if(finx) {
		fin.seekg(0,"end");
		len=fin.tellg();
		fin.seekg(0,"beg");
		fin>>header;
		try{if(!header.check(MyTop.Size())) throw string("Number of atoms in the pdb and trajectory files does not not match!");}
		catch(const string & s){cout << s<<endl;exit(-1);}
	}
	atm->initPDB(topPDB);
}
}
